# TODO: Add comment
# 
# Author: guochun
###############################################################################

source("sourceAll.R")
N=200
mf=function(mi,mj){
	return(exp(-abs(mi-mj)))
}
temp=numeric()
rep=100
marks=runif(N,1,10)
for(i in 1:rep){
	x=runif(N,0,10)
	y=runif(N,0,10)
	marks=sample(marks)
	X=ppp(x=x,y=y,window=owin(c(0,10),c(0,10)),marks=marks)
	re.k=Km(X,mf)
	temp=c(temp,re.k[,2])
}
dim(temp)=c(dim(re.k)[1],rep)
re.k[,2]=apply(temp,1,mean)
re.g=pcfConvert(re.k)
plot(x=re.g[,1],y=re.g[,2],type="l")
ep1=numeric()
for(i in 1:50000){
	tp=sample(marks,2)
	ep1[i]=mf(tp[1],tp[2])
}
abline(h=mean(ep1),col=2)
#looks as the same as expectation

#step 2: test the relationship with Ripley's K function
X2=rpoispp(function(x,y) { 6 * exp(-x/6)}, win=owin(c(0,10),c(0,10)))
marks(X2)=runif(X2$n,1,10)
mf2=function(mi,mj){
	return(rep(1,length(mj)))
}
re.k2=Km(X2,mf2)
re.k2[,2]=re.k2[,2]*mean(ep1)
re.g2=pcfConvert(re.k2)
temp=numeric()
rep=100
N=X2$n
marks=marks(X2)
for(i in 1:rep){
	X2$marks=sample(marks)
	re.k=Km(X2,mf)
	temp=c(temp,re.k[,2])
}
dim(temp)=c(dim(re.k)[1],rep)
re.k[,2]=apply(temp,1,mean)
re.g=pcfConvert(re.k)
plot(x=re.g[,1],y=re.g[,2],type="l")
lines(x=re.g[,1],y=re.g2[,2],col=2)
# also looks as expectation

#test the normal confidence interval
N=500
X=ppp(x=runif(N,0,10),y=runif(N,0,10),
		window=owin(c(0,10),c(0,10)),marks=runif(N,1,10))
re.k=Km(X,mf)
conf.simu=envelopeKm(X,mf,plot=TRUE)

reply.k=Km(X,mf2)
ep1=numeric()
marks=marks(X)
for(i in 1:50000){
	tp=sample(marks,2)
	ep1[i]=mf(tp[1],tp[2])
}
conf.theo.up=pcfConvert(data.frame(re.k[,1],re.k[,2]+1.96*reply.k[,2]*sd(ep1)/sqrt(40)))
conf.theo.low=pcfConvert(data.frame(re.k[,1],re.k[,2]-1.96*reply.k[,2]*sd(ep1)/sqrt(40)))
lines(x=re.k[,1],y=conf.theo.up[,2],col=3)
lines(x=re.k[,1],y=conf.theo.low[,2],col=4)

#there are some large bias, so this might not be such an easy job.


#let's generate six point patterns
# and test the correlation bewteen marks
N=500
Xa=ppp(x=runif(N,0,100),y=runif(N,0,100),window=owin(c(0,100),c(0,100)),marks=runif(N,1,30))
Xc=Xb=Xa
marks(Xb)=Xb$x
#Xc=rmpp(Xc,20,NA,c(1,30))
Xc=rpp(Xc$n,c(100,100))[owin(c(0,100),c(0,100))]

f=function(x,y) { 1/(x+1)^(0.5) }
N=N*5.7
Xd=ppp(x=runif(N,0,100),y=runif(N,0,100),window=owin(c(0,100),c(0,100)),marks=runif(N,1,30))
Xd=Xe=Xf=rthin(Xd,f)
marks(Xd)=runif(Xd$n,1,30)
marks(Xe)=Xe$x
Xf=rpp(N,c(100,100),TRUE,f)[owin(c(0,100),c(0,100))]

par(mfrow=c(2,3),mai=c(0.2,0.2,0.4,0.2))
plot(Xa,main="a")
plot(Xb,main="b")
plot(Xc,main="c")
plot(Xd,main="d")
plot(Xe,main="e")
plot(Xf,main="f")

windows()
par(mfrow=c(2,3),mai=c(0.5,0.5,0.4,0.5))
envelopeKm(Xa,mf,plot=TRUE,main="a",r=seq(0,10,0.5))
envelopeKm(Xb,mf,plot=TRUE,main="b",r=seq(0,10,0.5))
envelopeKm(Xc,mf,plot=TRUE,main="c",r=seq(0,10,0.5))
envelopeKm(Xd,mf,plot=TRUE,main="d",r=seq(0,10,0.5))
envelopeKm(Xe,mf,plot=TRUE,main="e",r=seq(0,10,0.5))
envelopeKm(Xf,mf,plot=TRUE,main="f",r=seq(0,10,0.5))


#test difference between two marked point patterns by the second null model
par(mfrow=c(1,2))
envelopeKm(Xc,Xe,mf,plot=TRUE,r=seq(0,10,0.5))
envelopeKm(Xc,Xf,mf,plot=TRUE,r=seq(0,10,0.5))

#emprical test
#two kind of questions can be asked, one is how the mark function distributed in space?
#the other one is that is there any spatial correlation between marks?
data=getData("bci5.txt")
#subdata=data[which(data$gx<=100 & data$gy >= 400 & !is.na(data$dbh) &!is.na(data$gx) &!is.na(data$gy)),]
spab=table(data$sp)
selsp=names(spab)[spab>200 & spab<800]
pdf("./all species2.pdf")
for (i in 11:length(selsp)){
	sp1=selsp[i]
	subdata=data[data$sp==sp1 & !is.na(data$dbh) &!is.na(data$gx) &!is.na(data$gy),]
	pd=ppp(x=subdata$gx,y=subdata$gy,window=owin(c(0,1000),c(0,500)),marks=subdata$dbh)
	re.g=try(envelopeKm(pd,mf,plot=TRUE,main=sp1))
}
dev.off()

#########33
data0=getData("bci0.txt")
data5=getData("bci5.txt")
data=data0
subdata0=data[which(data$gx<=100 & data$gy >= 400 & !is.na(data$dbh) &!is.na(data$gx) &!is.na(data$gy)),]
data=data5
subdata5=data[which(data$gx<=100 & data$gy >= 400 & !is.na(data$dbh) &!is.na(data$gx) &!is.na(data$gy)),]

pd0=ppp(x=subdata0$gx,y=subdata0$gy,window=owin(c(0,100),c(400,500)),marks=subdata0$dbh)
pd5=ppp(x=subdata5$gx,y=subdata5$gy,window=owin(c(0,100),c(400,500)),marks=subdata5$dbh)

re=envelopeKm(pd5,pd0,mf,r=seq(0,25,0.5))
